// Filename: FFTFS.CPP

/************   Fast Fourier Transform   ********************\

1993 (c) CDI UI, The USSR Academy of Sciences, Samara
FFT function for normal precision (symmetric).

\************************************************************/

#include "fft.h"

int fftfs( Complex *x, int Size, int DirInv ) {

  if (Size==1) return NoError;
  if (Size > FFTFMAXSIZE || Size <= 0) return InvalidSize;

  int k = 1, Pow2 = 0; 
  while ( k<Size ) { k *= 2; Pow2++; }
  if (k!=Size || (DirInv!=1 && DirInv!=-1)) return InvalidParams;

  for (int i=1; i<Size-1; i+=2) { x[i].re = -x[i].re; x[i].im = -x[i].im; } 

  int j=0;
  for (i=0; i<Size-1; i++) {
    Complex t;
    if (i<j) { t=x[j]; x[j]=x[i]; x[i]=t; }
    k=Size/2;
    while (k<=j) { j-=k; k/=2; }
    j+=k;
  }

  for (int l=0, ll=1; l<Pow2; l++, ll*=2) {
    Complex u = {1, 0}, w; 
    w.re =  cos(DirInv * M_PI/ll);
    w.im = -sin(DirInv * M_PI/ll);
    for (j=0; j<ll; j++) {
      float tr, ti;
      for (i=j; i<Size; i+=2*ll) {
        k=i+ll;
        tr = x[k].re * u.re - x[k].im * u.im;
        ti = x[k].re * u.im + x[k].im * u.re;
        x[k].re = x[i].re - tr; x[k].im = x[i].im - ti;
        x[i].re += tr; x[i].im += ti;
      }
      tr = u.re * w.re - u.im * w.im;
      ti = u.re * w.im + u.im * w.re;
      u.re = tr;  u.im = ti;
    }
  }

  for (i=1; i<Size-1; i+=2) { x[i].re = -x[i].re; x[i].im = -x[i].im; } 

  return NoError;
}
